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Abstract 

We consider an extension of Kuramoto's model of coupled phase oscillators where 
oscillator pairs interact with different strengths. When the coupling coefficient of 
each pair can be separated into two different factors, each one associated to an 
oscillator, Kuramoto's theory for the transition to synchronization can be explicitly 
generalized, and the effects of coupling heterogeneity on synchronized states can be 
analytically studied. The two factors are respectively interpreted as the weight of the 
contribution of each oscillator to the mean field, and the coupling of each oscillator 
to that field. We explicitly analyze the effects of correlations between those weights 
and couplings, and show that synchronization can be completely inhibited when they 
are strongly anti-correlated. Numerical results validate the theory, but suggest that 
finite-size effect are relevant to the collective dynamics close to the synchronization 
transition, where oscillators become entrained in synchronized frequency clusters. 
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1 Introduction 



Macroscopic periodic oscillations are ubiquitous in living organisms, and span 
a broad range of time scales, from fractions of a second -such as in neural 
and heart tissues- to days or weeks -such as in circadian rhythms and men- 
strual cycles. The possibility that this kind of macroscopic dynamics is the 
collective manifestation of mutual synchronization of microscopic oscillations 
was advanced by N. Wiener in the 1940s [1], and was later elaborated by 
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A. T. Winfree [2,3]. The basic mechanism assumed to be at work in this self- 
organization phenomenon is that a few elements in a population of interacting 
oscillators may synchronize if they have similar frequencies and their coupling 
is strong enough. Under suitable conditions, other oscillators may in turn be 
entrained by this "nucleation center" and, eventually, contribute to form a 
macroscopic oscillating cluster. Collective oscillations would thus emerge as 
the result of coherent microscopic dynamics, through a process not unlike a 
condensation phase transition [2]. 

This conceptual idea was realized into an explicit model by Y. Kuramoto [4] , 
who considered the collective dynamics of an ensemble of globally coupled 
phase oscillators. The state of each oscillator is defined by its phase (j)i{t) £ 
[0, 27i). Their joint evolution is governed by the equations 

K ^ 



{i — 1, . . . , A") where uji is the natural frequency of oscillator i, and K is 
the strength of their global (all-to-all) couphng. In this system, synchronized 
states can be characterized by the distribution of effective frequencies uj[. The 
effective frequency of an oscillator is defined as the time average of its phase 
velocity: 

1 ^ 

uj[=lim^-J^idt. (2) 





It is found that, for sufficiently large K, there is a group of oscillators whose 
effective frequencies are identical. The number of these synchronized oscillators 
increases with K. Synchronization implies as well a certain degree of order in 
the distribution of the individual oscillator phases. This is revealed by the 
complex order parameter 

1 ^ 

z{t) = a{t) exp[i$(i)] = ^ E exp mt)] ■ (3) 



Kuramoto showed that in the limit N ^ oo, and under the assumption that 
synchronized oscillators form a single group with effective frequency cul — ^, 
the collective amplitude a -which turns out to be independent of time- van- 
ishes for coupling intensities below a certain critical value K^, and is different 
from zero for K > [4]. The threshold Kc of the synchronization transition 
is determined by the distribution of natural frequencies ujf ensembles with 
wider frequency distributions require stronger coupling to become synchro- 
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nized. As a byproduct of the formulation, the synchronization frequency is 
also obtained, and = $(0) + Qt. 

The quantity z{t) can be interpreted as a mean field generated by the oscillator 
ensemble. In fact, a{t) and define, respectively, the collective amplitude 
and phase of the ensemble. The evolution equation (1) for each oscillator can 
be rewritten as 

(j)i^u;i + Ka sin(^ - (4) 



The couphng of oscillator i with the ensemble is thus represented as the in- 
teraction with a single, equivalent oscillator of phase weighted by the am- 
plitude a. 

The synchronization transition results from the competing effect of two fac- 
tors: while coupling induces the emergence of ordered dynamics, the hetero- 
geneity in the distribution of natural frequencies favours incoherent behaviour 
[5]. The aim of the present paper is to incorporate two additional sources of 
heterogeneity to Kuramoto's model. First, we consider that each oscillator i 
is coupled to the mean field z = (jexp(i$) with its own coupling intensity 
ki. Second, we assume that the individual contribution of each oscillator to 
the mean field is weighted by a factor qi. In contrast with natural frequen- 
cies, which determine the individual dynamics in the absence of coupling, the 
two new attributes ki and affect the way in which each oscillator interacts 
with the ensemble. Therefore, they define a heterogeneous interaction pattern 
underlying the system. 

Heterogeneity is a widespread feature in multi-agent systems. Constitutional 
differences between the members of such systems imply variations in the in- 
dividual dynamics, in the response to external influxes, in the interaction 
with other elements, etc. Those differences may also have non-trivial effects 
on the collective dynamics of the ensemble [3,6]. In connection with a more 
realistic representation of both natural and artificial systems, it is therefore 
desirable to generalize standard models in order to encompass heterogeneity. 
The present extension of the phase-oscillator model has the additional interest 
of being analytically tractable through a rather straightforward generalization 
of Kuramoto's theory. In fact, in the limit of an infinitely large ensemble, the 
synchronization transition, as well as the collective dynamical properties of 
the synchronized state, can be fully characterized in the frame of such gener- 
alization. 

In the next section, we discuss how Kuramoto's theory is extended to in- 
clude the form of heterogeneous coupling advanced above. We establish self- 
consistency equations for the collective amplitude and the synchronization 
frequency, and find the frequency distribution of non-synchronized oscillators. 
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In Section 3, the results of Section 2 are applied to analyze synchronization 
properties in several specific cases, paying special attention to the efi^ect of 
correlations between the different attributes of individual oscillators. In par- 
ticular, we show that anti-correlation between ki and Qi can completely sup- 
press synchronized states. Numerical results for finite systems are presented in 
Section 4, which validate analytical predictions and, at the same time, show 
that finite-size effects play an important role in the collective dynamics close 
to the synchronization transition. In this region, partial synchronization in 
the form of frequency clustering -not predicted by the analytical approach- 
is observed. Finally, our main results are summarized and commented in the 
last section. 



2 Synchronization transition for heterogeneous coupling 



Consider the collective dynamics of a population of N phase oscillators, whose 
phases 4>i{t) evolve according to 

1 ^ 



i = 1, . . . ,N. The coefficients Wij > weight the interaction of each oscillator 
pair. We assume that these coefficients can be expressed in the form of a 
product, Wij — kiQj, with > and > for all i. Equations (5) can be 
written as 

uji + kiU sin($ - (pi) (6) 



where a{t) and are the modulus and phase of the complex number [cf. 
Eq. (3)] 



1 



z{t) = a{t) exp[i$(i)] = ^ E ?i exp[i0,(i)]. (7) 



Kuramoto's model, described in the Introduction, is recovered for ki = K and 
Qi = 1 for all i. The special case ki = qi = Si for all i was studied by H. 
Daido [7,8]. However, inspired by disordered spin systems, he admitted Sj to 
be positive or negative. 

Equation (6) shows that, as in Kuramoto's model, the evolution of each oscil- 
lator can be though of as resulting from its interaction with the mean field z, 
generated by the oscillator ensemble. The factor ki can be interpreted as the 
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coupling of oscillator i with the mean field. Meanwhile, qi weights the contri- 
bution of the same oscillator to z. In view of this, we require that Qi satisfies 
the normalization condition 

1 

1^E?^ = 1- (8) 

i=l 



This choice does not imply loosing generality: if iV ^ Z^j = Q 7^ 1, we just 
redefine Qi/Q ^ Qi and Qki ki, leaving the evolution equations invariant. 

We focus on the limit of an infinitely large ensemble, N ^ 00, and describe 
the system in terms of continuous distributions. Specifically, the distribution 
of individual natural frequencies uji, couplings ki, and weights Qi, is given by 
a function G{u!,k,q). This distribution is normalized to unity, 

00 CXj OC' 

J dcu j dk j dq G{u;, k, q) = 1, (9) 
-00 



and, according to Eq. (8), it must satisfy the additional condition 

00 cx) 00 

j dw J dkj dqq G{u;, k, q) = 1. (10) 



-00 



As for the distribution of individual phases, the quantity n(0, t; /c, q)d(j) dk dq 
represents the fraction of oscillators with coupling in {k, k + dk) and weight 
in (g, q + dq), whose phases lie in (0, + dcj)) at time t. In terms of the density 
n{(f),t;k,q), Eq. (7) reads 

00 00 2n 

a"exp(i<I>) = j dk j q dq J dc/) n{(f),t; k, q) exp{i(f)). (11) 
00 



As in Kuramoto's theory, we assume that the density n{(l),t;k,q) has two 
contributions, n = nu + fig. The first contribution represents non-synchronized 
oscillators, which are uniformly distributed in phase, so that their density 
Hu does not depend on 4>. From Eq. (11), it is clear that non-synchronized 
oscillators do not contribute to the mean field, because the integration over 
vanishes. 

The second contribution represents a group of synchronized oscillators, which 
are assumed to possess a constant collective frequency fl. Their density has the 
form ns{(j), t; k, q) = ns{(f)—ftt; k, q). Replacing this in Eq. (11), the macroscopic 
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phase turns out to be ^{t) — $(0) + Qt, and 



oo oo 2n 



(jexp[i$(0)] = J dk J q dq J dip ns{'ip;k,q)exp{iip). (12) 







Within this Ansatz, thus, a is independent of time. With an appropriate choice 
of the origin of phases, we can fix $(0) = 0. 

Introducing the macroscopic phase = Qt in Eq. (6), and defining the 

individual deviation from $ as ipi{t) = 4>i{t) — we get 

ipi — LUi — fl — kia sin ipi. (13) 



In this equation for ipi, the collective amplitude a and the synchronization 
frequency fl are not known. However, we do know that, in the present for- 
mulation, they are independent of time. Therefore, Eq. (13) can be formally 
solved for generic constant values of a and fl. Once the deviations ipi -and, 
thus, the distribution of phases- have been found, the unknowns a and Q 
are calculated from Eq. (11). This self-consistent calculation of the collective 
properties of the oscillator ensemble is the core of Kuramoto's theory. 

When IcOi — il| < /cjcr, Eq. (13) has a stable fixed point at 



ipi = arcsm — , (14) 

kjCr 



with — 7r/2 < ipi < 7r/2. This fixed point stands for the (time-independent) 
asymptotic phase deviation of a synchronized oscillator of natural frequency 
uji and coupling ki with respect to the macroscopic phase ^{t) = Qt. Asymp- 
totically, the effective frequency of synchronized oscillators, Eq. (2), is uj^ — Q. 
For each value of the coupling, Eq. (14) relates the asymptotic phase of a 
synchronized oscillator with its natural frequency. Consequently, this equa- 
tion can be used to link the asymptotic density ns{ip;k,q) of synchronized 
oscillators to the distribution G{uj, k, q) of natural frequencies, couplings and 
weights. Taking into account the relation 

Usiip; k, q) dip dk dq = G{u!, k, q) dcu dk dq, (15) 
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which holds for ~ti /2 < i/j < 7r/2, and noting that Eq. (14) imphes the 
differential relation ka cos ip dip dk dq — duj dk dq, we find 



ns{ip;k,q) = < 



kaG{fl + ka sinip, k, q) cosip for — 7r/2 < ip < 7t/2, 



(16) 







otherwise. 



Replacing this form of the density of synchronized oscillators in Eq. (12), and 
separating real and imaginary parts, yields the self-consistency equations 



oo oo 



7r/2 

a — ajkdkjqdq J dtp G{fl + ka sin ■0, k, q) cos^ ip, 

-it/2 



(17) 



and 



oo oo 



7r/2 

/ 

-7r/2 



— ajkdkjqdq j dip G{Q + kasim/j,k,q) cosi/jsim/j, 



(18) 



for the collective amplitude a and the synchronization frequency Q, in terms 
of the distribution G{u}, k, q). Once these quantities are obtained, the density 
of synchronized oscillators can be explicitly calculated using Eq. (16). We 
postpone the discussion of the solutions to Eqs. (17) and (18), which define 
the synchronization properties of the oscillator ensemble, to the next section. 

Oscillators for which |cjj — i7| > fcjO" do not reach a stationary deviation from 
the macroscopic phase $ and, therefore, do not synchronize. In this case, the 
solution to Eq. (13) implies that the phase of a non-synchronized oscillator 
varies with time as 



(19) 



where ^i{t) is a 27r-periodic function of t. The effective frequency cul is given 

by 



\ 



ha 
uji — VL 



(20) 



For each value of /cj, this equation links the effective frequency u}[ with the 
natural frequency uji of each non-synchronized oscillator. It makes it pos- 
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sible to calculate the distribution G'{uj' . k, q) of effective frequencies, cou- 
plings and weights, in terms of G{uJ,k,q). Taking into account the relation 
G'{u!', k, q)duj' dk dq = G{u!, k, q)du! dk dq, and using Eq. (20) to obtain the 
ratio dcu/duj', we get 



G'ioj', k, q) = 



\UJ 



G 



n)^ + k^a^,k,q 



(21) 



This distribution of effective frequencies for non-synchronized oscillators can 
be explicitly calculated once a and fl have been found. 



3 Analysis of synchronization properties 



The self-consistency equations (17) and (18) define the collective macroscopic 
properties of the synchronized oscillators in terms of the distribution G{uj, k, q) 
of individual frequencies, couplings, and weights. A full analysis of the solu- 
tions for an arbitrary form of G{u!,k,q), taking into account any possible 
dependence on its three variables, is out of reach. In the following subsections, 
consequently, we study in detail a few representative special cases. 

It is however possible to begin by pointing out a few generic properties of 
Eqs. (17) and (18). First of all, these equations -and, as a matter of fact, the 
full formulation presented in Sect. 2- reduces to Kuramoto's theory when all 
oscillators have the same coupling K and unitary weight, i.e. when 

G{u, k, q) = g{u)6{k - K)6{q - 1), (22) 



as expected. In this limit, the only relevant individual attribute is the natural 
frequency Ui, with distribution g{u). 

For any form of G{uj,k,q), a trivial solution to Eqs. (17) and (18) is o" = 0. 
According to Eq. (16), this corresponds to a fully non-synchronized collec- 
tive state. To study synchronization features we look for solutions with non- 
vanishing mean field, a ^ 0. Moreover, as is known to happen in Kuramoto's 
theory, if the distribution G{uj, k, q) is an even function with respect to the 
variable cu around a certain value cuq, i.e. 

G{u;o + uj,k,q) = G{uJo-u},k,q), (23) 



Eq. (18) is satisfied by taking fl = ujq. In other words, for a symmetric dis- 
tribution of natural frequencies, the synchronization frequency coincides with 
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the center of the distribution. Now, since in Eqs. (5) natural frequencies are 
defined up to an arbitrary additive constant (corresponding to a phase rotation 
at constant angular velocity) we can always choose lvq = 0. In this case, the 
synchronization frequency vanishes. Under the above conditions, the problem 
reduces to Eq. (17) in the form 

oo oo ""/S 

1 = j k dk j q dq j d^jj G{ka simp, k,q) cos^ tjj. (24) 

-7r/2 



The synchronization threshold, at which the non-trivial solution reaches its 
lowest value a — 0, is defined by the condition 

oo oo 

J k dk J qG{0,k,q) dq^ -. (25) 



In the following, for the sake of conciseness, we deal with this specific situation. 

Before passing to the consideration of Eqs. (24) and (25) for some special forms 
of G{u!, k, q), let us point out an important property regarding the dependence 
of our problem on the variable q. Suppose that the weights qi are not correlated 
with the natural frequencies cui or with the couplings ki, so that we have 

G{uj,k,q) = G,{uj,k)G2{q). (26) 



Then, the distribution of weights becomes irrelevant to the problem. Indeed, 
under such condition, the integration over q in Eqs. (24) and (25) -as well as 
in Eqs. (17) and (18)- can be trivially performed, taking into account that 
Eq. (8) requires 

oo 

J q G2{q) dq = 1. (27) 





From then on, the distribution of weights plays no role in defining the macro- 
scopic synchronization properties of the ensemble. In the density of synchro- 
nized oscillators, Eq. (16), 6^2(5) appears as a trivial factor modulating the 
shape of ns{ip;k,q). This weak role of the variable q in the solution to the 
present problem is explained by noting that, in an infinitely large ensemble, 
an average quantity such as the mean field z does not change if each term in 
the average is weighted by a random coefficient, as long as the mean value 
of these coefficients equals unity. In our formulation, this property can be 
traced back to Eq. (6), which governs the evolution of each oscillator. Once 
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the mean field z has been introduced, the individual weights qi disappear from 
the calculation. 

This situation changes drastically if, on the other hand, there is a correlation 
between the weight and the couphng or the natural frequency of each oscillator. 
In this case, since the asymptotic phase of each synchronized oscillator depends 
on its coupling and frequency as given by Eq. (14), a correlation emerges 
between the weight and the phase. As a consequence, the mean field depends 
on how q is distributed, and on its correlation with lo and k. In view of this 
remark, our analysis of non-trivial distributions of q will focus on the case 
where the attributes /cj, and of each individual oscillator are mutually 
correlated. 



3.1 Uniform weights 



We address first the case where the weights are equal for all oscillators, qi — I 
for all i. In this situation, k, q) = Gi{u, k)S{q — 1), and Eq. (24) becomes 

oo t/2 

l^Jkdk J d7jjGi{kasmij;,k)cos^ijj. (28) 

-7r/2 



In Ref. [9], we performed a preliminary analysis of Eq. (28), paying particular 
attention to the case where couplings and natural frequencies are in turn un- 
correlated, Gi{u!, k) = g{uj)h{k). There, it was found that the synchronization 
threshold is given by the condition 

{k)^ [kh{k)dk^^ (29) 



In Kuramoto's theory, where all couplings are equal, ki = K = {k) for all 
i, the critical value of the coupling at which synchronization switches on is 
Kc = 2/7Tg{0). Therefore, Eq. (29) establishes that, in the case where couplings 
are distributed (and uncorrelated to natural frequencies), synchronization be- 
comes possible when the average couphng reaches Kuramoto's threshold Kc- 
Close to the synchronization transition, the collective amplitude behaves as 

{{k)-K,y/^ (30) 



nK,{k-^)\g"{0)\ 



where g"{0) is the second derivative of the distribution of natural frequencies at 
the center of the distribution, a; = 0. As in Kuramoto's theory, this is assumed 
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to be a negative number, corresponding to a maximum in the distribution. 

More generally, if the joint distribution of natural frequencies and couplings 
behaves as Gi{uj,k) — Go{k) — 7(^)0;^ close to a; = 0, the synchronization 
transition takes place when 

[kGo{k)dk^-. (31) 

J TT 




Close to the transition, the solution of Eq. (28) is approximately given by 



^ kGo{k)dk - 1 



(32) 



The possibility of having synchronization is determined both by the number 
of oscillators with natural frequencies at the center of the distribution and by 
their coupling distribution. In fact, Eq. (31) can be rewritten as 

2 

{k)Q go^ -, (33) 

TT 



where {k)o is the average coupling of the oscillators with cu — 0, and go — 
fo° Go{k)dk is the density of these oscillators. To become synchronized, a 
low number of oscillators at the center of the distribution requires that their 
average coupling is high, and vice versa. 



3.2 Uniform couplings 



When the couplings of all oscillators are equal, ki = K for all i, the distribution 
of natural frequencies, couplings, and weights is factorized as G{u!,k,q) = 
G2{u;, q)S{k - K), and Eq. (24) reads 

00 7r/2 

l = Kjqdq J dip G2{Kasmip,q)cos^ij. (34) 

-7r/2 



As discussed above, unless the function G2{uj,q) establishes a correlation be- 
tween weights and natural frequencies, the problem would be equivalent to 
Kuramoto's case. 



11 



Even if there is a correlation between weights and natural frequencies, this 
case of uniform couplings has the same mathematical form as Kuramoto's 
problem. In fact, Eq. (34) can be rewritten as 

7r/2 

1 = K j dip g{Kasinip)coB^il), (35) 

— 7r/2 



with 



g{uj) = j q G2{uJ,q) dq. 



(36) 



Equation (35) coincides with the equation for a of an ensemble with uniform 
couplings and weights, and with an effective frequency distribution g{uj). 

To evaluate the effect of correlations between weights and natural frequencies, 
we consider the extreme case where the weight qi of each oscillator is given 
by a function of its frequency cjj, say, qi = 6{uJi). In this case, G2(uj,q) = 
g{uj)6[q — 6{uj)]. Equation (10) requires g{uj)6{uj)duj = 1. The effective 
distribution of frequencies in Eq. (35) turns out to be g{uj) = 9{uj) g{uj), so 
that the condition for the synchronization transition reads 



Tim ^0{0)g{0)- 



This result shows that, with respect to the case of uncorrelated weights, the 
threshold couphng for synchronization decreases when the weights of the oscil- 
lators at the center of the frequency distribution are larger than the average, 
^(0) > 1, i.e. when the contribution of those oscillators to the mean field is 
relatively strong. On the contrary, if the weight of those oscillators in relatively 
small, a stronger coupling is necessary to induce synchronization. 



3.3 Weight-coupling correlation 



A more interesting situation arises when both weights and couplings adopt 
non- uniform values over the ensemble. In this shown below, the cor- 

relation between ki and qi can have drastic effects on the collective behaviour 
of oscillators, to the point of suppressing synchronization even in the presence 
of arbitrarily strong couplings. 

To simplify the discussion, we assume that individual natural frequencies are 
not correlated with weights or couplings, so that G{LJ,k,q) = g{Lj)Gs{k,q). 
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As for the weight-coupling correlation, wc consider the extreme case in which 
Qi is a given function of ki, qi — Q{ki), so that G^{k,q) — h{k)S[q — Q{k)]. 
Condition (10) requires that 

oo 

J h{k)e{k)dk = 1, (38) 





and the threshold for synchronization is determined by the equation 

f k h(k) e(k) dk^ (39) 



Suppose now that the product h{k)Q{k) satisfies Eq. (38), and compare Eq. 
(39) with Eq. (29), which determines the synchronization threshold for the 
case of uniform weights. U q = Q{k) is an increasing function of k -i.e., 
if weights and couplings are positively correlated- and for a given value of 
g{0), the synchronization condition is expected to hold in the present case for 
relatively lower couphngs. On the other hand, if Q{k) decreases with k -i.e., 
if weights and couplings are anti-correlated- relatively larger couplings will 
be required to reach the synchronization threshold. It may even happen, as 
we show below, that a sufficiently strong anti-correlation between weights and 
couplings suppresses the occurrence of synchronization even when arbitrarily 
large couplings are present in the ensemble. 

From the viewpoint of the dynamical roles of weights and couplings, this effect 
of positive or negative correlation between ki and qi is interpreted as follows. 
The emergence of synchronized behaviour is facilitated if the oscillators whose 
coupling with the mean field z is stronger arc in turn those whose contribution 
to z is relatively larger. This kind of feedback between oscillators with large 
ki and qi enhances the development of coherent states. On the other hand, if 
the mean field is dominated by oscillators whose coupling with the ensemble 
is weak, synchronization can be inhibited or even fully suppressed. 

To illustrate these phenomena, we explicitly work out a case with a particularly 
simple distribution of couplings, which is later used again in our numerical 
simulations of Section 4. Let us consider a uniform distribution, h{k) = 
for < A; < -fCmax, and h{k) — otherwise. As for the weights, we take 

A K 

q = Q{k) = ^ exp(Afc), (40) 

exp(AKmax) - 1 

where the pre- factor has been chosen to satisfy Eq. (38). The weight-coupling 
correlation is positive for A > 0, because q grows with k. For A < 0, on the 
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other hand, q and k are anti-correlated. In the hmit A = 0, the case of uniform 
weights is reobtained. 



With these choices, Eq. (39) yields 
1 - (1 - Ai^max) exp{XK 



A[exp(Ai^^ax) - 1] 7r^(0)' 



(41) 



which can be interpreted as an equation for the maximum couphng -ft'max at 
the synchronization transition. For smaU A, the solution to this equation reads 

K ~ ^ - (42) 



As expected, with respect to the case of uniform weights -where the synchro- 
nization threshold occurs at i^max = 4/77(7(0)- a lower value of -ft'max is required 
if weights and couplings are positively correlated (A > 0). On the other hand, 
if A < and the correlation is negative, couplings must be larger to induce 
synchronization. 

Note now that, for negative A, the left-hand side of Eq. (41) is an increasing 
function of i^max, which asymptotically approaches the value |A|^^ as -ft'max 
grows to infinity. This imphes that, if A < Ac = — 7rg'(0)/2, there will be no 
finite value of -ft'max such that Eq. (41) is satisfied. As A approaches the criti- 
cal value Ac from above, in fact, the solution for -ft'max diverges to infinity. In 
other words, if the anti-correlation between weights and couplings is strong 
enough, even arbitrarily large couplings are unable to induce the synchro- 
nization transition. For A < Ac and any K^^x the only possible value for the 
collective amplitude is the trivial one, a — 0. 



4 Numerical results. Frequency clustering 



In this section, we present a series of results from the numerical resolution of 
Eqs. (5). The aim of these numerical calculations is two-fold. On the one hand, 
they validate our main analytical results over wide regions of the relevant pa- 
rameter spaces. On the other, they show that a systematic departure from the 
analytical results is observed, in finite systems, close to the synchronization 
transition. As advanced in Ref. [9], this departure is due to the occurrence 
of a regime of frequency clustering, presumably ascribable to finite-size fiuc- 
tuations, and disregarded in the analytical approach. Part of our numerical 
calculations is thus devoted to illustrate and characterize this regime, where 
a portion of the ensemble segregates into several groups of mutually synchro- 
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nized oscillators. Within each group, all oscillators share the same effective 
frequency. 

We solve Eqs. (5) by means of a standard Euler algorithm, for ensembles of 
N — 10^ phase oscillators. In order to emphasize differences with the case 
of homogeneous global coupling, we focus on variations in the distribution of 
couplings ki and weights qi, and on their possible correlation, while we assume 
that natural frequencies are uncorrelated to coupling and weights. The distri- 
bution of natural frequencies is a Gaussian, g{u!) = exp(— a;^/2)/\/27r. Results 
are typically obtained from averages over 100 realizations of the ensemble, with 
couplings, weights, and natural frequencies drawn anew from their respective 
distributions. 

A crucial point in our numerical calculations is the identification of mutu- 
ally synchronized oscillators or, cquivalently, of synchronized clusters. This is 
achieved by comparing the distribution of natural frequencies cui and effective 
frequencies u'^ [9]. In a given realization of the ensemble, natural frequencies 
are drawn at random from the above-referred Gaussian distribution, and effec- 
tive frequencies are calculated from a numerical approximation to the integral 
in Eq. (2). The distributions of cUj and uj'^ are constructed as histograms with 
columns of width A = 10~^. In the histogram of natural frequencies we identify 
the highest column, whose height we denote by h. In our calculations, typical 
values of h are 10 to 15. A column in the histogram of effective frequencies is 
considered to belong to a synchronized cluster if its height is above h. A syn- 
chronized cluster is defined as a set of contiguous columns higher than h whose 
nearest columns to the left and to the right are lower than h. Synchronized 
oscillators are those oscillators belonging to synchronized clusters. 



4-1 Uncorrelated couplings and weights 

First, we have measured the total number Ng of synchronized oscillators for 
two uncorrelated distributions of couplings and weights. The upper panel of 
Fig. 1 shows the fraction of synchronized oscillators, = Ng/N, in the case of 
constant couplings, ki = K for all i, as a function of K. Solid dots correspond 
to the case of constant weights, Qi = 1 for all i. This case coincides with 
Kuramoto's model, which we use as a reference system. Open dots, on the 
other hand, stand for ensembles where the weights are drawn at random from 
a uniform distribution over the interval (0, 1), and then normalized to satisfy 
Eq. (8). 

The full curve represents the analytical estimate for the fraction of synchro- 
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K 



Fig. 1. Upper panel: Numerical results for the fraction Ug of synchronized oscillators 
in the case of constant couplings, ki = K for all z, with constant weights, = 1 for 
all i (solid dots), and with uniformly distributed weights (open dots). Lower panel: 

The same, for the case of uniformly distributed couplings, < fcj < Xmax- Full and 
dotted curves stand for the analytical calculation of Ug and the collective amplitude 
(7, respectively. 



nized oscillators, defined from our analytical results as 



27r 



dk j dq j dij) ns{^\k,q), 



(43) 



while the dotted curve is the analytical result for the collective amplitude 
(J, obtained from Eq. (24). Both Us and cr are different from zero above Ku- 
ramoto's critical coupling Kc — 2/7rg'(0) = 2J2/t: 1.596. 
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The agreement between numerical and analytical results is very good for 
K > Kc- Below the synchronization transition, on the other hand, we find 
a noticeable difi^erence which, as advance above, is associated with the forma- 
tion of synchronized clusters. This phenomenon is illustrated later. Note also 
that, in agreement with our discussion on the irrelevance of the distribution 
of weights when they are uncorrelated to couplings, solid and open dots above 
the transition fall over the same curve. 

The lower panel of Fig. 1 displays analogous results for the case where cou- 
plings are uniformly distributed in the interval (0,i^^max)- The corresponding 
distribution is h{k) — K~l^ for < A; < i^^max, and h{k) — otherwise, and 
the synchronization transition takes place at -ft'max = ~ 3.192. Again, 
the agreement between numerical and analytical results is very good above 
the transition, and the same kind of deviation is observed below it. Also, as 
expected, the distribution of weights qi has no effect on the fraction of syn- 
chronized oscillators. 

A suitable illustration of the state of the system close to the synchronization 
transition, where analytical and numerical results differ from each other, is 
provided by a plot of the effective frequencies cj' versus the natural frequencies 
uji of all oscillators in a given realization of the ensemble. In this kind of 
plot, a horizontal array of dots -corresponding to oscillators with different 
natural frequencies but the same effective frequency- reveals the presence of 
a synchronized cluster. Figure 2 displays such plots for the four cases already 
considered in Fig. 1. In the upper line, couplings are constant, ki = K = 1.5 
for all i, while in the lower line they are uniformly distributed between and 
-f^max = 3. In both cases, thus, the system is just below the synchronization 
threshold. The appearance of synchronized clusters near the center of the 
frequency distribution is apparent in the four plots. 

Figure 2 shows again that, for a given distribution of couplings, there is lit- 
tle qualitative difference between the cases where the weights qi are either 
constant or distributed. On the other hand, we find that for oscillators with a 
given natural frequency, the effective frequencies are essentially identical in the 
case of constant couplings, but exhibit a substantial dispersion for distributed 
couplings. This difference can be understood by extrapolating our analytical 
results for the synchronization regime to the present situation, just below the 
transition. In fact, Eq. (20) implies that, for a given natural frequency, the 
effective frequency depends on ki, so that we expect to have the same value 
of uj^ in the case of uniform coupling, and different values if couplings are 
distributed. It also predicts that the effective frequency of a given oscillator 
should be smaller than its natural frequency when the latter is larger than 
the synchronization frequency Q, and vice versa. We see from Fig. 2 that our 
numerical calculations, for which Q 0, agree with such prediction. 
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Fig. 2. The effective frequency as a function of the natural frequency uji of indi- 
vidual oscillators in single realizations for constant coupling, ki = 1.5 for all i, with 
constant weights (upper right), and uniformly distributed weights (upper left), and 
for uniformly distributed couplings, < ki < 3, with constant weights (lower right), 
and uniformly distributed weights (lower left). Only the central region of the fre- 
quency distributions is displayed. 



Our previous results for the case of constant weights [9], which we do not 
reproduce here, suggest that the regime of frequency clustering is a finite- 
size effect due to fluctuations in the distribution of natural frequencies. An 
unusual accumulation of oscillators in a given frequency interval can trigger 
local entrainment for couplings below the transition to synchronization. Nu- 
merical results show that clustering becomes less important, with a smaller 
fraction of entrained oscillators, as the ensemble grows in size. In the limit of 
an infinitely large system, fluctuations should average out and synchronization 
should switch on at the analytically predicted threshold, with the formation of 
a macroscopic cluster. However, the fact the frequency clustering is still con- 
spicuous in a relatively large ensemble as in our calculations seems to point 
out that this regime plays a relevant role in the collective dynamics of finite 
systems. 
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Fig. 3. Synchronization region (gray shaded), obtained analyticahy for anti-cor- 
relatcd couphngs and weights, Eq. (40), in the (i^maxj l^^l) parameter plane. The 
inserts show plots of individual effective frequencies versus natural frequencies (cf. 
Fig. 2), for single realizations of the system at the indicated points of the parameter 
plane. Scales in all the inserts vary from —0.4 to 0.4 on both axes. 

4-2 Weight- coupling anti- correlation 



We turn now the attention to the case discussed in Section 3.3, where couplings 
and weights are correlated. As above, couplings are uniformly distributed on 
(0, ii'max), and the weight of each oscillator is given as function of its coupUng 
through Eq. (40). We focus on the more interesting case of weight-coupling 
anti-correlation, where synchronization suppression is possible, so that we take 
A < 0. Results are thus presented in terms of the positive parameter |A|. We 
recall that, in our analysis of infinitely large ensembles, Eq. (41) predicts the 
relation between A and Kmax at the synchronization threshold. Figure 3 shows 
the region of the (i^^max, |A|) parameter plane where synchronized states are 
possible. Synchronization is suppressed for small couplings and strong weight- 
coupling anti-correlation. In particular, it cannot occur for any |A| if i^max < 
2/ng{0) = 4.^2/^ ^ 3.192, or for any K^ax if |A| > ng{0)/2 = ^ 0.627. 

Inserts in Fig. 3 illustrate the relation of effective and natural frequencies at 
several points in the parameter plane. 

For small values of |A|, the main effect of weight- coupling anti-correlation is 
to shift the synchronization transition to higher couplings. This is shown in 
Fig. 4, where the fraction of synchronized oscillators is plotted as a function of 
Kuiax for four values of |A|. Full lines correspond to the analytical prediction. 
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Fig. 4. The fraction of synchronized oscillators for the case of anti-correlated 
couplings and weights, as a function of Kmax, for four values of |A|. Full curves 
stand for the analytical prediction, and dotted curves are plotted as a guide to the 
eye in the region of frequency clustering. 

which again gives a very good description of the numerical results above the 
transition, at least for |A| = 0.2 and 0.4. Dotted lines are plotted as a guide 
to the eye in the region preceding the synchronization transition. For |A| = 
0.6, the systems is practically at the hmit where synchronization becomes 
inhibited for any -fTmax- The analytical transition point is strongly shifted 
to the right, but a considerable fraction of the ensemble is synchronized in 
clusters well below the transition. Note that a discrepancy between numerical 
and analytical results persists above the critical point. Finally, for |A| = 1, 
synchronization should be completely suppressed for infinitely large systems. 
In our numerical calculations, a small fraction of the ensemble -which grows 
with Xmax" is however synchronized. 

It turns out that the quantities that characterize synchronization in the clus- 
tering regime -such as the fraction of synchronized oscillators, and the number 
of clusters- depend strongly on the value of i^max- This is illustrated, for in- 
stance, by the two inserts of Fig. 3 corresponding to parameters just outside 
the synchronization region. For i^max = 4, clusters are small and relatively 
sparse, and most effective frequencies are close to the corresponding natural 
frequencies (all dots are near the diagonal). For Xmax = 10, on the other hand, 
a substantial fraction of the ensemble around the center of the frequency dis- 
tribution is entrained into clusters, and form massive groups extended over 
considerable intervals of natural frequencies. This effect can be understood 
taking into account that a fluctuation in the distribution of frequencies is 
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Fig. 5. Upper panel: Fraction Ug of synchronized oscillators as a function of |A|, for 
Kmax = 4. The full curve is the analytical prediction, and the dotted curve has been 
added as a guide to the eye in the clustering regime. Lower panel: The corresponding 
number of clusters, C. The dotted curve is a spline approximation. 

more likely to give rise to a synchronized cluster if the involved oscillators 
have, on the average, larger couplings. 

A more quantitative illustration of this dependence on -ft^max is provided by the 
numerical results presented in Figs. 5 and 6. The upper panel of Fig. 5 shows 
the fraction Us of synchronized oscillators as a function of |A|, for K^aax — 4. 
The curve corresponds to the analytical prediction. For |A| larger than the 
critical value, Ug decreases rapidly, such that only 2% of the ensemble remains 
clustered for |A| ~ 1. In the lower panel, we plot the corresponding number 
of clusters, C. Below the transition, in the synchronization region, C ~ 1. 
Around the transition point, the number of clusters increases abruptly, and 
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reaches a maximum just below C — 20, before beginning to decline for larger 
|A|. 



For Xmax = 10, the behaviour of both rig and C is qualitatively similar, but 
important quantitative differences are apparent. First, the fraction of syn- 
chronized oscillators in the clustering region is much higher and its decay is 
slower, with approximately 10% of the ensemble still entrained into clusters for 
|A| ~ 1. Second, the number of clusters at the transition grows above C = 60, 
and the subsequent decay is extremely slow. For values of |A| twice as large as 
the critical point, C remains practically as large as just above the transition. 



5 Conclusion 

In this paper, we have considered an ensemble of coupled phase oscillators 
where the strength of the interaction is generally different for each oscillator 
pair. This heterogeneity is represented by symmetric positive interaction coef- 
ficients Wij and, together with the distribution of natural frequencies uji, adds 
diversity to the ensemble. Our main goal has been to show that Kuramoto's 
theory for the synchronization transition of an infinitely large ensemble of 
globally coupled phase oscillators can be extended to the case of heteroge- 
neous coupling, yielding exact results when the interaction coefficients can be 
factorized as Wij = kiqj. The factors ki and qi turn out to be, respectively, 
the coupling strength of oscillator i with the mean field generated by the en- 
semble, and the weight with which the same oscillator contributes to that 
mean field. The three attributes that characterize each oscillator, cjj. ki. and 
Qi, are distributed over the ensemble according to a function G{u}, k, q) which, 
in general, introduces statistical correlations between them. 

In the analysis of our results, we have paid particular attention to the effects 
of correlations between natural frequencies cUj, couplings fcj, and weights q^. 
As it may have been expected, if the couplings and/or the weights of the 
oscillators which trigger synchronization -around the center of the frequency 
distribution- are relatively small, an overall larger coupling will be necessary 
to effectively have synchronized states. The effect of correlations between cou- 
plings and weights, on the other hand, is less obvious. We have found that 
a sufficiently strong anti-correlation between ki and g,j is able to completely 
inhibit synchronization, in the sense that synchronized states are suppressed 
even in the presence of arbitrarily large couplings. In other words, if the mean 
field is strongly dominated by oscillators whose coupling is weak, and vice 
versa, the ensemble may be unable to display coherent collective oscillations. 

As a validation of the analytical formulation, we have performed numerical 
calculations for ensembles of 10^ oscillators, with several combinations of the 
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Fig. 6. As in Fig. 5, for K^^^ = 10. 

distribution of couplings and weights. Numerical and analytical results are in 
good overall agreement, except for a systematic departure in the parameter re- 
gions where the synchronization transition takes place. In these zones, numer- 
ical results reveal the existence of a regime of partial synchronization, where 
oscillators become entrained into several internally synchronized clusters. As 
the transition is approached, these clusters grow in size and, at the same 
time, they progressively aggregate with each other. Eventually, they collapse 
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into a single cluster, which can be identified with the macroscopic fraction 
of synchronized oscillators predicted by the analytical formulation. Previous 
numerical results on a special case of the present system [9] , seem to indicate 
that frequency clustering is a finite-size effect which disappears for infinitely 
large systems -where, precisely, the analytical approach is formulated. In finite 
ensembles, clusters may emerge from localized fluctuations in the distribution 
of frequencies and couplings, which trigger the mutual synchronization of a 
few oscillators. 

The regime of clustering preceding the transition to synchronization found 
in our system, is reminiscent of a similar phenomenon well-known to occur 
in ensembles of coupled chaotic elements [10,11,12,5]. In the present case, 
heterogeneity in the ensemble replaces chaotic dynamics as the factor which 
counteracts the effects of coupling and gives origin to the transition. In spite 
of the fact that it may disappear in the thermodynamical limit, clustering 
remains the richest dynamical regime of oscillator ensembles because of its 
complexity and diversity [13]. 

In connection with chaotic systems, it would be interesting to study whether 
the introduction of coupling heterogeneity in ensembles of chaotic dynamical 
elements have effects similar to those described here for the synchronization 
transition of phase oscillators. More complex individual dynamics are also ex- 
pected to give rise to new collective phenomena, such as amplitude effects. 
Distributed couplings ki and weights Qi can be straightforwardly incorporated 
to standard models of interacting chaotic systems [5]. For instance, linear het- 
erogeneous coupling between chaotic elements whose individual dynamics 
is governed by the equation x = f(x) may be introduced as 

Xi = f(xO + A;i(X-Xi), (44) 

i = 1, . . . , N, with X = A^~^ J2i Qi^i- Finally. let us point out that Kuramoto's 
theory has been extended, at least partially, to study the synchronization prop- 
erties of ensembles formed by other kinds of dynamical elements, in particular, 
by active rotators [14,15]. Consideration of such ensembles with the kind of 
heterogeneous coupling studied here would be a natural continuation of the 
present work. 
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